library(signal)
library(cba)
library(reshape2)
library(XML)

#******************************
# G L O B A L   V A R S
#******************************


#******************************
# C O N S T A N T S
#******************************

SAMPLE_SIZE=10000
sensors=list(inside=c('10A9F57B28040','10C717C280AF','10D147C28094'),outside=c('1095E87B28040'))
file = 'c:\\igloo.log'

#******************************
# F U N C T I O N S
#******************************

queryNOAA <- function() {
	node=parseXMLAndAdd(readLines(url('http://www.weather.gov/xml/current_obs/KHEF.xml'))[-1])
	for (obs in c('observation_time','latitude','longitude',
		'station_id','temp_c','pressure_in','relative_humidity','dewpoint_c')
	) {
		cat(sprintf("%s = %s\n",obs,getXMLNodeText(node,obs)))
	}
}

getXMLNodeText <- function(node,name) {
	xmlValue(getNodeSet(node,sprintf("/tmp/current_observation/%s",name,sep=''))[[1]])
}

stopWatch <- local({
	watch.time = NA
	function() {
		if (is.na(watch.time)) {
			watch.time <<- Sys.time()
			cat("Starting watch at",format(watch.time,format="%m/%d %H:%M"))
		} else {
			new.time <- Sys.time()
			cat(sprintf("Stopwatch delta %.2f s", new.time - watch.time))
			watch.time <<- new.time
		}
	}
})

plotTemps <- function(data,title) {
	ylim=c(
		max(-10,min(data,na.rm=T)),
		min(35,max(data,na.rm=T))
	)
	matplot(data,typ='l',ylim=ylim,lty=1,xaxt='n',ylab='Temp (C)',main=title)
	drawAxes(stamps.lt,format="%b %d %Hh",temp.mult=10,hour.mult=6,ylim=ylim)
	#legend(x=0,y=20,legend=unlist(sensors),cex=0.5,col=1:3,lty=1,lwd=1)
}

#http://knol.google.com/k/fast-fourier-transform-and-power-spectrum-using-r-software
plotFFT <- function(data, title) {
	# calculate fft of data (this gives the complex transform)
	fourier = fft(data)
	 
	# extract the power spectrum (power is sometimes referred to as "magnitude")
	magnitude = sqrt(Re(fourier)*Re(fourier)+Im(fourier)*Im(fourier))
	 
	# extract the phase 
	phase = atan(Im(fourier)/Re(fourier))
	
	# select only first half of vectors
	magnitude_firsthalf <- magnitude[1:(length(magnitude)/2)]
	phase_firsthalf<-phase[1:(length(magnitude)/2)]
	 
	# generate x-axis
	x.axis <- 1:length(magnitude_firsthalf)/length(magnitude)
	 
	# plot the power spectrum
	plot(x=x.axis,y=magnitude_firsthalf,type='l',log='xy')
}

drawAxes <- function(s,format="%m/%d %H:%M",temp.mult=10,hour.mult=2,ylim) {
	# vline indexes
	hour.threshold.idxs=( which( (s$hour - c(s$hour[-1],NA) !=0)) + 1 )
	hour.grid.idxs=hour.threshold.idxs[which((1:length(hour.threshold.idxs) %% hour.mult) == 1)] #--- every other hour.mult hours
	labels=format(s[hour.grid.idxs],format=format)
	axis(1, at=hour.grid.idxs, labels=labels,cex.axis=0.7,las=3)
	temp.lines=(ylim[1]:ylim[2])[which(ylim[1]:ylim[2] %% temp.mult == 0)]
	abline(h=temp.lines,col='gray')
	abline(v=hour.grid.idxs,col='gray')
	day.lines=which((s$yday - c(s$yday[-1],NA))<0)
	abline(v=day.lines,col='blue')
} 

getRawData <- function() {
	stopWatch()
	cat("Reading igloo buffer file")
	arduino.all = readLines(file)[which(count.fields(file,sep='|')==7)]
	stopWatch()
	return(arduino.all)
}
	
cleanData <- function(arduino.all,sample.size) {
	arduino.sensor = arduino.all[grep('\\|DS18S20\\|',arduino.all,value=F)]
	arduino = sample(arduino.sensor,size=sample.size,replace=F)
	arduino = sort(arduino)
	
	cat("Parsing")
	
	d=read.table(text=arduino,sep='|',row.names=NULL,header=F,as.is=T,fill=T,
		colClasses=c('POSIXct','character','character','character','character','character','numeric')
		)
		
	names(d)=c('ts','master','port','slave','sensor','model','temp')
	#d$ts=as.POSIXct(round(d$ts))
	
	invalid_row_ids=which(!d$sensor %in% unlist(sensors))
	if (length(invalid_row_ids)>0) {
		cat(sprintf("Discarding %i rows from invalid sensors\n",length(invalid_row_ids)))
		invalid_row_ids
		d=d[-invalid_row_ids,]
	}
	
	summary(d)
	return(d)
}

meltData <- function(d) {
	stopWatch()
	cat("Melting")
	dm=(melt(d,id.vars=c('ts','sensor'),measure.vars=c('temp')))
	dm.summary=summary(dm$value)
	dm=dcast(dm,ts~sensor)
	
	summary(dm)
	
	i=1
	for (sensor in unique(d$sensor)) {
		if (i==1) {
			w = which(!is.na(dm[,sensor]))
			out = dm[w,c('ts',sensor)]
		} else {
			w = which(!is.na(dm[,sensor]))[1:nrow(out)]
			out = cbind(out,dm[w,sensor])
			names(out)[ncol(out)]=sensor
		}
		i = i + 1
	}
	return(out)
}

#******************************
# M A I N
#******************************

igloo <- function () {
	stopWatch()
	data.all=getRawData()
	data.clean=cleanData(data.all,sample.size=SAMPLE_SIZE)
	data=meltData(data.clean)
	summary(data)
	stopWatch()
	
	cat("Converting dates")
	stamps.lt=as.POSIXlt(round(data$ts),origin='1970-01-01') # drop second fractions to align to second resolution
	
	# ***************************************
	# ******* G R A P H I C S ***************
	# ***************************************
	
	graphics.off()
	par(mfrow=c(2,2))
	
	title=sprintf("Raw signal (%i/%i samples)",SAMPLE_SIZE,length(data.all))
	plotTemps(data[,unlist(sensors)],title=title)
	
	filter.size=30
	title=sprintf("Linear f: %i (%.2f min)",filter.size,difftime(data[filter.size,'ts'],data[1,'ts'],units='mins'))
	outf=stats::filter(data[,unlist(sensors)], filter=rep(1/filter.size,filter.size), method = c("convolution"),sides = 1, circular = FALSE)
	plotTemps(outf,title=title)
	
	title=sprintf("Rolling median (k=7)")
	out.rmedian=rollapply(data=data[,unlist(sensors)],width=7,FUN=median)
	plotTemps(out.rmedian,title=title)
	
	title=sprintf("Rolling mean (k=10)")
	out.rmean=rollmean(x=data[,unlist(sensors)],k=10)
	plotTemps(out.rmean,title=title)
	
	#plotFFT(data[,3],'FFT')
	#time.axis=data.frame(ts=as.POSIXlt(seq(min(stamps.lt),max(stamps.lt),by="sec")))
	#data.proj=merge(x=time.axis,y=data,all.x=T)
	#sparse time series too large at second grain -- up to minute?
	#space between readings seems x50 (100, 150), is it? outliers?
	#messed up dcast (melt) when ts was made ct instead of lt?
	
	#library(zoo)
}
